Mouse muscle stem cell engraftment: Bulk RNA-seq gene differential analysis

Génome Québec Innovation Centre, McGill University, Montreal, Canada

Prepare analysis workflow

# knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# options(
#   readr.show_progress = FALSE,
#   digits = 2,
#   scipen = 8,
#   future.globals.maxSize = +Inf
# )
knitr::opts_knit$set(root.dir = getwd())

library(tidyverse)
library(tximport)
library(DESeq2)
library(vsn)
library(cowplot)
library(glmpca)
library(apeglm)
library(rhdf5)
library(EnhancedVolcano)
library(pheatmap)
library(ensembldb)
library(scater)
library(PoiClaClu)
library(RColorBrewer)
require(edgeR)
theme_set(theme_bw())

Auxiliary functions

get_filtered_genes_index <- function(scaled_tpm, meta_dt, min_count = 100, min_total = 150) {
  design_matrix_groups <-
    model.matrix(
      ~ condition:batch,
      meta_dt %>%
        column_to_rownames("sample")
    )

  keep <- edgeR::filterByExpr(as.matrix(scaled_tpm[, -c(1, 2)]),
    design_matrix_groups,
    min.count = min_count,
    min.total.count = min_total
  )
  return(keep)
}


get_deseq_data <- function(target_tissue, design_matrix, gene_metadata, meta_dt, data_dir, tx2knownGene) {
  samples_dir <- file.path(
    data_dir,
    meta_dt %>%
      pull(name)
  )
  cdna_dir <- file.path(samples_dir)
  cdna_files <- file.path(cdna_dir, "abundance.h5")
  names(cdna_files) <-
    meta_dt %>%
    pull(sample)

  # bias corrected counts without an offset.  Corrects for potential differential isoform usage

  txi_scaledTPM <- tximport(cdna_files,
    type = "kallisto",
    tx2gene = tx2knownGene,
    countsFromAbundance = c("scaledTPM")
  )


  scaled_tpm <-
    txi_scaledTPM$counts %>%
    as_tibble(rownames = "gene_id")


  scaled_tpm <-
    left_join(scaled_tpm,
      gene_metadata %>%
        dplyr::select(gene_name, gene_id),
      by = "gene_id"
    ) %>%
    relocate(gene_name, .after = gene_id)

  write_tsv(scaled_tpm, paste0(target_tissue, "_engraftment_scaled_tpm.tsv"))

  # for deseq2
  txi <- tximport(cdna_files,
    type = "kallisto",
    tx2gene = tx2knownGene,
    countsFromAbundance = c("no")
  )


  dds <-
    DESeqDataSetFromTximport(
      txi,
      meta_dt %>%
        column_to_rownames("sample"),
      design = design_matrix
    )

  keep <- get_filtered_genes_index(scaled_tpm, meta_dt)
  table(keep)
  dds <- dds[keep, ]
  design_matrix_reduced <-
    model.matrix(
      ~1,
      meta_dt %>%
        column_to_rownames("sample")
    )
  colnames(design_matrix_reduced) <- c("intercept")
  dds <- DESeq(dds, test = "LRT", reduced = design_matrix_reduced)

  counts_normalized <-
    counts(dds, normalized = T) %>%
    as_tibble(rownames = "gene_id") %>%
    left_join(.,
      gene_metadata %>%
        dplyr::select(gene_name, gene_id),
      by = "gene_id"
    ) %>%
    relocate(gene_name, .after = gene_id)

  write_tsv(counts_normalized, paste0(target_tissue, "_engraftment_counts_filtered_normalized.tsv"))

  return(dds)
}

# We use apeglm shrinkage to preserve the size of large LFC and compute s-values (the estimated rate of false sign). Shrinkage is useful for ranking and visualization, without the need for arbitrary filters on low count genes. https://doi.org/10.1093/bioinformatics/bty895

# The local false sign rate FSR is defined as the posterior probability that the posterior mode (MAP) of beta (log2FC) is of a biologically significant effect size (i.e. abs(beta) > lfcThreshold )

get_lfc <- function(dds, coef_name, lfc_threshold, svalue_threshold) {
  lfc <-
    lfcShrink(dds,
      coef = coef_name,
      type = "apeglm",
      svalue = T,
      lfcThreshold = lfc_threshold
    )

  plotMA(lfc, ylim = c(-22, 8))
  abline(h = c(-1, 1), col = "dodgerblue", lwd = 2)

  lfc <-
    lfc %>%
    as_tibble() %>%
    dplyr::select(-baseMean) %>%
    mutate(
      keep = ((abs(log2FoldChange) > lfc_threshold) & (svalue < svalue_threshold))
    ) %>%
    set_names(c(
      paste0(
        c("lfc_", "lfc_se_", "svalue_", "keep_"),
        coef_name
      )
    ))

  return(lfc)
}

get_lfc_estimates <- function(dds, gene_metadata, coef_names, lfc_threshold, svalue_threshold) {
  
  lfcs <-
    map_dfc(coef_names,
      get_lfc,
      dds = dds,
      lfc_threshold = lfc_threshold,
      svalue_threshold = svalue_threshold
    )

  #coef_names_batch <- c(paste0(c("intercept", coef_names), "_batch"), paste0("SE_", c("intercept", coef_names), "_batch"))

  lfcs <-
    bind_cols(
      mcols(dds)[, c("baseMean", "baseVar", "dispOutlier", "maxCooks", "intercept", "SE_intercept", coef_names)] %>%
        as_tibble(rownames = "gene_id") %>%
        dplyr::select(gene_id, everything()),
      lfcs
    )



  lfcs <-
    lfcs %>%
    dplyr::select(
      gene_id, baseMean,
      intercept, paste0("lfc_", coef_names),
      paste0("svalue_", coef_names), everything()
    )

  lfcs <-
    lfcs %>%
    left_join(., gene_metadata, by = "gene_id") %>%
    dplyr::select(gene_name, everything()) %>%
    mutate(
      gene_name =
        scater::uniquifyFeatureNames(
          ID = gene_id,
          names = gene_name
        )
    )
  return(lfcs)
}


plot_volcano <- function(data, xvar, yvar, title, lfc_threshold, svalue_thresh, figures_dir) {
  gv <-
    EnhancedVolcano(data,
      lab = data %>% pull(gene_name),
      xlab = bquote(~ italic(Moderated) ~ Log[2] ~ FC),
      ylab = bquote(~ -Log[10] ~ italic(svalue)),
      x = xvar,
      y = yvar,
      col = c("grey30", "forestgreen", "red2", "royalblue"),
      pCutoff = svalue_thresh,
      FCcutoff = lfc_threshold,
      #  legend = c("NS", "Log2FC", "svalue", "svalue & Log2FC"),
      legendLabels = c(
        "NS",
        expression(Log[2] ~ FC),
        "svalue",
        expression(s - value ~ and ~ log[2] ~ FC)
      )
    )

  ggsave(file.path(figures_dir, paste0(title, "_volcano.pdf")), gv)
  return(gv)
}

plot_heatmap <- function(target_genes, title, meta_col, data, show_genes = F, figures_dir) {
  df <-
    left_join(target_genes %>%
      dplyr::select(gene_id),
    data,
    by = "gene_id"
    ) %>%
    dplyr::select(-gene_id) %>%
    column_to_rownames("gene_name") %>%
    as.matrix()

  gh <-
    pheatmap(df,
      color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 11, name = "Spectral")))(100),
      scale="row",
      cluster_rows = T,
      show_rownames = show_genes,
      treeheight_row = 0,
      fontsize_col = 4,
      fontsize_row = 6,
      angle_col = "45",
      cluster_cols = FALSE,
      annotation_col = meta_col
    )
  ggsave(file.path(figures_dir, paste0(title, "_heatmap_expression.pdf")), gh)

  return(gh)
}

plot_pca <- function(vsd, title, figures_dir) {
  pcaData <-
    plotPCA(vsd,
      intgroup = c("time", "age", "tissue", "batch", "donor"),
      returnData = TRUE
    )
  percentVar <- round(100 * attr(pcaData, "percentVar"))

  gg_plot <-
    ggplot(
      pcaData %>%
        mutate(condition = paste(age, time, sep = "_")),
      aes(x = PC1, y = PC2, color = time, shape = age, label = name)
    ) +
    geom_point(size = 4) +
    geom_line(aes(group = donor), colour = "grey") +
    #  coord_fixed() +
    geom_text(check_overlap = T, angle = 0, size = 3, vjust = 3, nudge_y = 1.5) +
    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
    coord_fixed() +
    ggtitle("PCA with VST data")
  ggsave(file.path(figures_dir, paste0(title, "_pca_bulk_rnaseq.pdf")), gg_plot)
  gg_plot
}

plot_gene <- function(gene, data, figures_dir) {
  df_plot <-
    data %>%
    dplyr::filter(gene_id == gene)

  gene_name <- df_plot$gene_name[1]
  p <- ggplot(
    df_plot,
    aes(x = time, y = lognorm_counts)
  ) +
    geom_point(aes(color = donor)) +
    geom_line(aes(group = donor)) +
    labs(title = gene_name) +
    facet_grid(batch ~ age)
  ggsave(file.path(figures_dir, paste0(gene_name, "_expression_lognorm_counts.pdf")), p)
  return(p)
}

Set filepaths

data_dir <- "./data/kallisto_out"
gene_metadata_filepath <- "./data/metadata/mm10_ens97_gene_meta.txt"
metadata_filepath <- "./data/metadata/samples_metadata_leg.tsv"

Get tx id versions

library(AnnotationHub)

Loading required package: BiocFileCache
Loading required package: dbplyr

Attaching package: 'dbplyr'

The following objects are masked from 'package:dplyr':

    ident, sql


Attaching package: 'AnnotationHub'

The following object is masked from 'package:Biobase':

    cache

ah <- AnnotationHub()

DEPRECATION: As of AnnotationHub (>2.23.2), default caching location has changed.
  Problematic cache: /home/rick/.cache/AnnotationHub
  See https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/TroubleshootingTheCache.html#default-caching-location-update
snapshotDate(): 2021-05-18

edb <- ah[["AH73905"]]

loading from cache

seqlevelsStyle(edb) <- "UCSC"
# standard_chroms <- standardChromosomes(edb, species = organism(edb))
tx_meta <-
  transcripts(edb,
    columns = c(
      "tx_id_version",
      "gene_id",
      "gene_name",
      "gene_id_version"
    ),
    return.type = "DataFrame"
  )
tx2knownGene <-
  tx_meta %>%
  as_tibble() %>%
  dplyr::select(tx_id_version, gene_id)

meta_dt <-
  read_tsv(file.path(metadata_filepath)) %>%
  # arrange(tissue, age, time) %>%
  mutate(
    condition = paste(age, time, sep = "_"),
    age = factor(age, levels = c("young", "aged")),
    time = factor(time, levels = c("T0", "T21")),
    batch = factor(batch, levels = c(1, 2))
  )

Rows: 23 Columns: 8── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): name, sample, time, tissue, age, donor, type
dbl (1): batch
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

meta_dt

gene_metadata <-
  read_csv(gene_metadata_filepath) %>%
  dplyr::select(-gene_id_version)

Rows: 55487 Columns: 8── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): gene_name, gene_id, seqnames, gene_biotype, gene_id_version, description
dbl (2): gene_len, perc_gene_gc
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

gene_metadata <-
  left_join(tx_meta %>% as_tibble() %>%
    group_by(gene_id) %>%
    dplyr::count(gene_name) %>%
    dplyr::select(-n),
  gene_metadata %>%
    dplyr::select(-gene_name),
  by = "gene_id"
  ) %>%
  mutate(
    gene_name =
      scater::uniquifyFeatureNames(
        ID = gene_id,
        names = gene_name
      )
  )
gene_metadata

design_matrix <-
  model.matrix(
    ~ time * age * batch,
    meta_dt %>%
      column_to_rownames("sample")
  )
colnames(design_matrix) <- c(
  "intercept",
  "engraf",
  "age",
  "intercept_batch",
  "niche",
  "engraf_batch",
  "age_batch",
  "niche_batch"
)

design_matrix <- design_matrix[, c(1, 2, 3, 5, 4, 6, 7, 8)]
design_matrix

            intercept engraf age niche intercept_batch engraf_batch age_batch niche_batch
yng_t0_01           1      0   0     0               0            0         0           0
yng_t0_02           1      0   0     0               0            0         0           0
yng_t0_03           1      0   0     0               0            0         0           0
yng_t0_04           1      0   0     0               1            0         0           0
yng_t0_05           1      0   0     0               1            0         0           0
yng_t0_06           1      0   0     0               1            0         0           0
yng_t21_01          1      1   0     0               0            0         0           0
yng_t21_02          1      1   0     0               0            0         0           0
yng_t21_04          1      1   0     0               1            1         0           0
yng_t21_05          1      1   0     0               1            1         0           0
yng_t21_06          1      1   0     0               1            1         0           0
aged_t0_01          1      0   1     0               0            0         0           0
aged_t0_02          1      0   1     0               0            0         0           0
aged_t0_03          1      0   1     0               1            0         1           0
aged_t0_04          1      0   1     0               1            0         1           0
aged_t0_05          1      0   1     0               1            0         1           0
aged_t0_06          1      0   1     0               1            0         1           0
aged_t21_01         1      1   1     1               0            0         0           0
aged_t21_02         1      1   1     1               0            0         0           0
aged_t21_03         1      1   1     1               1            1         1           1
aged_t21_04         1      1   1     1               1            1         1           1
aged_t21_05         1      1   1     1               1            1         1           1
aged_t21_06         1      1   1     1               1            1         1           1

Load data

dds <- get_deseq_data("leg", design_matrix, gene_metadata, meta_dt, data_dir, tx2knownGene)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
summarizing abundance
summarizing counts
summarizing length
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
summarizing abundance
summarizing counts
summarizing length
using counts and average transcript lengths from tximport
using supplied model matrix
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
2 rows did not converge in beta, labelled in mcols(object)$fullBetaConv. Use larger maxit argument with nbinomLRT

dds

class: DESeqDataSet 
dim: 17153 23 
metadata(1): version
assays(6): counts avgTxLength ... H cooks
rownames(17153): ENSMUSG00000000001 ENSMUSG00000000028 ... ENSMUSG00000118487 ENSMUSG00000118488
rowData names(33): baseMean baseVar ... deviance maxCooks
colnames(23): yng_t0_01 yng_t0_02 ... aged_t21_05 aged_t21_06
colData names(8): name time ... batch condition

plotDispEsts(dds)

res <- results(dds)
W <- res$stat
maxCooks <- apply(assays(dds)[["cooks"]], 1, max)
idx <- !is.na(W)
plot(rank(W[idx]), maxCooks[idx],
  xlab = "rank of Wald statistic",
  ylab = "maximum Cook's distance per gene",
  ylim = c(0, 5), cex = .4, col = rgb(0, 0, 0, .3)
)
m <- ncol(dds)
p <- 3
abline(h = qf(.99, p, m - p))

plot(res$baseMean + 1, -log10(res$pvalue),
  log = "x", xlab = "mean of normalized counts",
  ylab = expression(-log[10](pvalue)),
  ylim = c(0, 30),
  cex = .4, col = rgb(0, 0, 0, .3)
)

lfc_threshold <- 1
svalue_thresh <- 0.05
coef_names <- colnames(design_matrix)[-1]
lfcs <-
  get_lfc_estimates(
    dds,
    gene_metadata,
    coef_names,
    lfc_threshold,
    svalue_thresh
  )

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

lfcs

    ggplot(lfcs, aes(x = lfc_niche_batch, y = lfc_age_batch)) +
    geom_point(size = 1)

  lfcs <-
    lfcs %>%
    mutate(
      total_batch = abs(intercept_batch) + abs(niche_batch) + abs(age_batch) + abs(engraf_batch),
      keep_niche_b = (keep_niche & !keep_niche_batch),
      keep_age_b = (keep_age & !keep_age_batch),
      keep_engraf_b = (keep_engraf & !keep_engraf_batch),
      EAN = paste0(
        as.integer(keep_engraf_b),
        as.integer(keep_age_b),
        as.integer(keep_niche_b)
      ),
    )

table(lfcs$EAN)

  000   001   010   011   100   101   110   111 
15374   608   314   118   351    78   240    70 

write_tsv(lfcs, "leg_two_batches_moderated_lfc_1_svalue_05_marked.txt")

summary_dt <-
  colSums(assay(dds)) %>%
  tibble::enframe(name = "sample", value = "counts")
summary_dt

scaled tpm

scaled_tpm <- read_tsv("leg_engraftment_scaled_tpm.tsv")

Rows: 36558 Columns: 25── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr  (2): gene_id, gene_name
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, aged_t2...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

counts_normalized <- read_tsv("leg_engraftment_counts_filtered_normalized.tsv")

Rows: 17153 Columns: 25── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr  (2): gene_id, gene_name
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, aged_t2...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

counts_normalized

counts_normalized_long <-
  counts_normalized %>%
  pivot_longer(cols = yng_t0_01:aged_t21_06, names_to = "sample", values_to = "counts") %>%
  left_join(., meta_dt %>% dplyr::select(sample, time, age, donor, batch, condition), by = "sample") %>%
  mutate(lognorm_counts = log1p(counts))

lfcs_hits <-
  lfcs %>%
  dplyr::filter(lfcs$EAN != "000")
lfcs_hits

targets <- c("Ccnd1", "Sparc",  "Col3a1", "Mt2", "Ccl11", "Col4a2")
 lfcs_genes <-
  lfcs %>%
  dplyr::filter(gene_name %in% targets)
  lfcs_genes

lfcs_genes_dt<-
  lfcs_genes %>%
  dplyr::select(gene_name, EAN, intercept:lfc_niche, intercept_batch:niche_batch, svalue_engraf:svalue_niche, keep_engraf, keep_age, keep_niche,  total_batch)
lfcs_genes_dt

figures_dir <- "./figures/hits"
map(lfcs_genes %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

Hits by group

target_genes_leg_001 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("001")) %>%
  arrange(svalue_niche)
target_genes_leg_001

figures_dir <- "./figures/hits/leg_001"
pp1 <- map(target_genes_leg_001 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_011 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("011")) %>%
  arrange(-baseMean)
target_genes_leg_011

figures_dir <- "./figures/hits/leg_011"
pp2 <- map(target_genes_leg_011 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_101 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("101"))
target_genes_leg_101

figures_dir <- "./figures/hits/leg_101"
pp3 <- map(target_genes_leg_101 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_111 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("111"))
target_genes_leg_111

figures_dir <- "./figures/hits/leg_111"
pp4 <- map(target_genes_leg_111 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_100 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("100"))
target_genes_leg_100

figures_dir <- "./figures/hits/leg_100"
pp5 <- map(target_genes_leg_100 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_010 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("010"))
target_genes_leg_010

figures_dir <- "./figures/hits/leg_010"
pp6 <- map(target_genes_leg_010 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_110 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("110"))
target_genes_leg_110

figures_dir <- "./figures/hits/leg_110"
pp7 <- map(target_genes_leg_110 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)

Saving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche", "svalue_niche", "niche", lfc_threshold, svalue_thresh, figures_dir)

One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_age", "svalue_age", "age", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_engraf", "svalue_engraf", "engraf", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche_batch", "svalue_niche_batch", "niche_batch", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_age_batch", "svalue_age_batch", "age_batch", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_engraf_batch", "svalue_engraf_batch", "engraf_batch", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_intercept_batch", "svalue_intercept_batch", "intercept_batch", lfc_threshold, svalue_thresh, figures_dir)

One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] AnnotationHub_3.0.2         BiocFileCache_2.0.0         dbplyr_2.1.1                edgeR_3.34.1                limma_3.48.3                RColorBrewer_1.1-2          PoiClaClu_1.0.2.1           scater_1.20.1              
 [9] scuttle_1.2.1               SingleCellExperiment_1.14.1 ensembldb_2.16.4            AnnotationFilter_1.16.0     GenomicFeatures_1.44.2      AnnotationDbi_1.54.1        pheatmap_1.0.12             EnhancedVolcano_1.10.0     
[17] ggrepel_0.9.1               rhdf5_2.36.0                apeglm_1.14.0               glmpca_0.2.0                cowplot_1.1.1               vsn_3.60.0                  DESeq2_1.32.0               SummarizedExperiment_1.22.0
[25] Biobase_2.52.0              MatrixGenerics_1.4.3        matrixStats_0.61.0          GenomicRanges_1.44.0        GenomeInfoDb_1.28.4         IRanges_2.26.0              S4Vectors_0.30.2            BiocGenerics_0.38.0        
[33] tximport_1.20.0             forcats_0.5.1               stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_2.1.2                 tidyr_1.1.4                 tibble_3.1.6               
[41] ggplot2_3.3.5               tidyverse_1.3.1            

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                    tidyselect_1.1.1              RSQLite_2.2.9                 grid_4.1.2                    BiocParallel_1.26.2           munsell_0.5.0                 ScaledMatrix_1.0.0           
  [8] ragg_1.2.1                    preprocessCore_1.54.0         withr_2.4.3                   colorspace_2.0-2              filelock_1.0.2                ggalt_0.4.0                   knitr_1.37                   
 [15] rstudioapi_0.13               Rttf2pt1_1.3.9                labeling_0.4.2                bbmle_1.0.24                  GenomeInfoDbData_1.2.6        farver_2.1.0                  bit64_4.0.5                  
 [22] coda_0.19-4                   vctrs_0.3.8                   generics_0.1.1                xfun_0.29                     R6_2.5.1                      ggbeeswarm_0.6.0              rsvd_1.0.5                   
 [29] locfit_1.5-9.4                bitops_1.0-7                  rhdf5filters_1.4.0            cachem_1.0.6                  DelayedArray_0.18.0           assertthat_0.2.1              vroom_1.5.7                  
 [36] promises_1.2.0.1              BiocIO_1.2.0                  scales_1.1.1                  beeswarm_0.4.0                gtable_0.3.0                  beachmat_2.8.1                ash_1.0-15                   
 [43] affy_1.70.0                   rlang_1.0.0                   genefilter_1.74.1             systemfonts_1.0.3             splines_4.1.2                 lazyeval_0.2.2                rtracklayer_1.52.1           
 [50] extrafontdb_1.0               broom_0.7.12                  BiocManager_1.30.16           yaml_2.2.2                    modelr_0.1.8                  backports_1.4.1               httpuv_1.6.5                 
 [57] extrafont_0.17                tools_4.1.2                   affyio_1.62.0                 ellipsis_0.3.2                jquerylib_0.1.4               Rcpp_1.0.8                    plyr_1.8.6                   
 [64] sparseMatrixStats_1.4.2       progress_1.2.2                zlibbioc_1.38.0               RCurl_1.98-1.5                prettyunits_1.1.1             viridis_0.6.2                 haven_2.4.3                  
 [71] fs_1.5.2                      magrittr_2.0.2                reprex_2.0.1                  mvtnorm_1.1-3                 ProtGenerics_1.24.0           mime_0.12                     hms_1.1.1                    
 [78] evaluate_0.14                 xtable_1.8-4                  XML_3.99-0.8                  emdbook_1.3.12                readxl_1.3.1                  gridExtra_2.3                 compiler_4.1.2               
 [85] biomaRt_2.48.3                bdsmatrix_1.3-4               maps_3.4.0                    KernSmooth_2.23-20            crayon_1.4.2                  htmltools_0.5.2               later_1.3.0                  
 [92] tzdb_0.2.0                    geneplotter_1.70.0            lubridate_1.8.0               DBI_1.1.2                     proj4_1.0-11                  MASS_7.3-55                   rappdirs_0.3.3               
 [99] Matrix_1.4-0                  cli_3.1.1                     pkgconfig_2.0.3               GenomicAlignments_1.28.0      numDeriv_2016.8-1.1           xml2_1.3.3                    annotate_1.70.0              
[106] vipor_0.4.5                   bslib_0.3.1                   XVector_0.32.0                rvest_1.0.2                   digest_0.6.29                 Biostrings_2.60.2             rmarkdown_2.11               
[113] cellranger_1.1.0              DelayedMatrixStats_1.14.3     restfulr_0.0.13               curl_4.3.2                    shiny_1.7.1                   Rsamtools_2.8.0               rjson_0.2.21                 
[120] lifecycle_1.0.1               jsonlite_1.7.3                Rhdf5lib_1.14.2               BiocNeighbors_1.10.0          viridisLite_0.4.0             fansi_1.0.2                   pillar_1.6.5                 
[127] lattice_0.20-45               ggrastr_1.0.1                 KEGGREST_1.32.0               fastmap_1.1.0                 httr_1.4.2                    survival_3.2-13               interactiveDisplayBase_1.30.0
[134] glue_1.6.1                    png_0.1-7                     BiocVersion_3.13.1            bit_4.0.4                     stringi_1.7.6                 sass_0.4.0                    blob_1.2.2                   
[141] textshaping_0.3.6             BiocSingular_1.8.1            memoise_2.0.1                 irlba_2.3.5